/*
 * topological-hex, a program to compute combinatorial hexahedral meshes.
 *
 * Copyright (C) <2018> <Université catholique de Louvain (UCL), Belgique>
 *
 * List of the contributors to the development, description and complete
 * License: see LICENSE file.
 *
 * This program (topological-hex) is free software:
 * you can redistribute it and/or modify it under the terms
 * of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program (see COPYING file).  If not,
 * see <http://www.gnu.org/licenses/>.
 */

#include <stdlib.h>
#include <string.h>

#include "symmetry.h"

static void invert_permutation(vertex_index *inv,
                               const vertex_index *permutation,
                               vertex_index n);

static
void init_incremental_symmetries(solver *solver);

error_code solver_set_symmetries(solver *solver,
                                 size_t num_symmetries,
                                 uint32_t num_boundary_vertices,
                                 const vertex_index *symmetries) {
  ALIGNED vertex_index *group;

  if (posix_memalign((void**)&group, ALIGNMENT,
                     sizeof(*solver->boundary_group) *
                     num_boundary_vertices * num_symmetries) != 0)
    return OUT_OF_MEMORY;

  if (solver->boundary_group)
    free(solver->boundary_group);

  solver->boundary_group = group;
  solver->num_symmetries = num_symmetries;
  solver->num_boundary_vertices = num_boundary_vertices;

  memcpy(solver->boundary_group, symmetries,
         sizeof(*solver->boundary_group) *
         num_boundary_vertices * num_symmetries);

  if (solver->flags & SOLVER_SHELLABLE_ONLY) {
    ALIGNED uint8_t *symmetry_data;
    if (posix_memalign((void**)&symmetry_data, ALIGNMENT,
                       sizeof(*symmetry_data) *
                       solver_incremental_symmetry_data_size(
                         solver->num_vertices, solver->num_hexes,
                         solver->num_symmetries)) != 0)
      return OUT_OF_MEMORY;

    if (solver->symmetry_data)
      free(solver->symmetry_data);

    solver->symmetry_data = symmetry_data;
    init_incremental_symmetries(solver);
  }

  return SUCCESS;
}

static
void init_incremental_symmetries(solver *solver) {
  size_t offset = sizeof(uint64_t);

  *(uint64_t*)solver->symmetry_data = 0;

  for (size_t i = 0; i < solver->num_symmetries; i++) {
    *(uint64_t*)solver->symmetry_data |= (1ul << i);

    memcpy(solver->symmetry_data + offset,
           solver->boundary_group + i * solver->num_boundary_vertices,
           sizeof(*solver->symmetry_data) * solver->num_boundary_vertices);
    for (vertex_index j = solver->num_boundary_vertices; j < solver->num_vertices; j++)
      solver->symmetry_data[offset + j] = VERTEX_SET_MAX;
    offset += solver->num_vertices;

    invert_permutation(solver->symmetry_data + offset,
                       solver->boundary_group + i * solver->num_boundary_vertices,
                       solver->num_boundary_vertices);
    for (vertex_index j = solver->num_boundary_vertices; j < solver->num_vertices; j++)
      solver->symmetry_data[offset + j] = VERTEX_SET_MAX;
    offset += solver->num_vertices;

    offset = (offset + alignof(uint64_t) - 1) & ~((size_t)(alignof(uint64_t) - 1));

    size_t n = sizeof(uint64_t) * ((solver->num_hexes + 63) / 64) + sizeof(uint64_t);
    memset(solver->symmetry_data + offset, 0, n);

    offset += n;
  }
}

static quad_index find_unvisited_quad(vertex_index *quads, uint32_t num_quads,
                                      bool *seen,
                                      vertex_index *permutation,
                                      vertex_index *reordered_quad) {
  for (quad_index i = 0; i < num_quads; i++) {
    if (seen[i]) continue;
    for (quad_vertex_index j = 0; j < 4; j++) {
      if (permutation[quads[4 * i + j]] != VERTEX_SET_MAX &&
          permutation[quads[4 * i + (j + 1) % 4]] != VERTEX_SET_MAX) {
        for (quad_vertex_index k = 0; k < 4; k++)
          reordered_quad[k] = quads[4 * i + (j + k) % 4];
        return i;
      }
    }
  }

  return NO_QUAD;
}

static quad_index find_quad_with_edge(vertex_index *quads, uint32_t num_quads,
                                      bool *seen,
                                      vertex_index *edge,
                                      vertex_index *reordered_quad) {
  for (quad_index i = 0; i < num_quads; i++) {
    if (seen[i]) continue;

    for (quad_vertex_index j = 0; j < 4; j++) {
      if (quads[4 * i + j] == edge[0] &&
          quads[4 * i + (j + 1) % 4] == edge[1]) {
        for (quad_vertex_index k = 0; k < 4; k++)
          reordered_quad[k] = quads[4 * i + (j + k) % 4];
        return i;
      }
      else if (quads[4 * i + j] == edge[1] &&
               quads[4 * i + (j + 1) % 4] == edge[0]) {
        reordered_quad[0] = quads[4 * i + (j + 1) % 4];
        reordered_quad[1] = quads[4 * i + j];
        reordered_quad[2] = quads[4 * i + (j + 3) % 4];
        reordered_quad[3] = quads[4 * i + (j + 2) % 4];
        return i;
      }
    }
  }

  return NO_QUAD;
}

bool find_symmetry(vertex_index *quads, uint32_t num_quads,
                   uint32_t num_vertices,
                   quad_index q_a, quad_index q_b,
                   uint8_t offset, bool flip,
                   bool *seen_a, bool *seen_b,
                   vertex_index *permutation,
                   vertex_index *inv_permutation) {
  for (vertex_index i = 0; i < num_vertices; i++)
    permutation[i] = inv_permutation[i] = VERTEX_SET_MAX;

  for (quad_index i = 0; i < num_quads; i++)
    seen_a[i] = seen_b[i] = false;

  quad_vertex_index flip_perm[] = {0, 3, 2, 1};

  for (quad_vertex_index i = 0; i < 4; i++) {
    vertex_index in_a, in_b;

    in_a = quads[4 * q_a + i];

    if (flip)
      in_b = quads[4 * q_b + flip_perm[(i + offset) % 4]];
    else
      in_b = quads[4 * q_b + ((i + offset) % 4)];

    permutation[in_b] = in_a;
    inv_permutation[in_a] = in_b;
  }

  seen_a[q_a] = seen_b[q_b] = true;

  vertex_index v_a[4], v_b[4];
  while ((q_a = find_unvisited_quad(quads, num_quads,
                                    seen_a, inv_permutation, v_a)) !=
         NO_QUAD) {
    vertex_index edge[] = {inv_permutation[v_a[0]], inv_permutation[v_a[1]]};

    q_b = find_quad_with_edge(quads, num_quads, seen_b, edge, v_b);
    if (q_b == NO_QUAD)
      return false;

    seen_a[q_a] = seen_b[q_b] = true;

    for (quad_vertex_index i = 2; i < 4; i++) {
      if (!extend_permutation(permutation, inv_permutation, v_b[i], v_a[i]))
        return false;
    }
  }

  return true;
}

error_code solver_compute_symmetries(solver *solver) {
  bool seen_a[solver->max_quad_id];
  bool seen_b[solver->max_quad_id];

  vertex_index permutation[solver->max_vertex_id];
  vertex_index inv_permutation[solver->max_vertex_id];

  vertex_index *permutations = malloc(sizeof(permutation) * 16);
  size_t capacity = 16;
  size_t size = 0;

  if (!permutations)
    return OUT_OF_MEMORY;

  for (quad_index i = 0; i < solver->max_quad_id; i++) {
    for (quad_vertex_index offset = 0; offset < 4; offset++) {
      for (int flip = 0; flip < 2; flip++) {
        if ((solver->flags & SOLVER_SHELLABLE_ONLY) == 0) {
          if (i == 0 && flip == 0 && offset == 0) continue;
          if (flip == 0) continue;
        }

        if (find_symmetry(solver->quads,
                          solver->max_quad_id, solver->max_vertex_id,
                          0, i, offset, flip, seen_a, seen_b,
                          permutation, inv_permutation)) {
          if (size == capacity) {
            size_t new_capa = (capacity * 3) / 2;
            vertex_index *buffer = realloc(permutations, sizeof(permutation) *
                                           new_capa);
            if (!buffer) {
              free(permutations);
              return OUT_OF_MEMORY;
            }

            permutations = buffer;
            capacity = new_capa;
          }

          memcpy(permutations + size * solver->max_vertex_id,
                 permutation, sizeof(permutation));
          size++;
        }
      }
    }
  }

  error_code status = solver_set_symmetries(
    solver, size, solver->max_vertex_id, permutations);
  free(permutations);

  return status;
}

static error_code no_goods_reserve(solver *solver, size_t to_add) {
  if (solver->no_goods_size + to_add > solver->no_goods_capacity) {
    size_t new_capa = (3 * solver->no_goods_capacity) / 2;
    if (solver->no_goods_size + to_add > new_capa)
      new_capa = solver->no_goods_size + to_add;

    void *new_buffer;
    if (posix_memalign(&new_buffer, ALIGNMENT,
                       new_capa * sizeof(*solver->no_goods)) != 0)
      return OUT_OF_MEMORY;

    if (solver->no_goods) {
      memcpy(new_buffer, solver->no_goods,
             solver->no_goods_size * sizeof(*solver->no_goods));
      free(solver->no_goods);
    }

    solver->no_goods = new_buffer;
    solver->no_goods_capacity = new_capa;
  }

  return SUCCESS;
}

error_code solver_start_branch(solver *solver, no_goods_delta *delta) {
  hex_index hex_id = solver->hex_id;

  if (solver->num_symmetries == 0)
    return SUCCESS;

  size_t old_size = solver->no_goods_size;

  if (solver->flags & SOLVER_SHELLABLE_ONLY) {
    if (solver->hex_id != 0) {
      if (no_goods_reserve(solver, 9) != SUCCESS)
        goto out_of_memory;

      solver->no_goods[solver->no_goods_size++] = NO_GOODS_HEX;

      memcpy(solver->no_goods + solver->no_goods_size,
             solver->hexes + 8 * (solver->hex_id - 1),
             8 * sizeof(*solver->hexes));

      solver->no_goods_size += 8;
    }

    if (no_goods_reserve(solver, 1) != SUCCESS)
      goto out_of_memory;

    solver->no_goods[solver->no_goods_size++] = NO_GOODS_DOWN;

    delta->size = solver->no_goods_size - old_size;

    return SUCCESS;
  }

  if (solver->no_goods_size == 0) {
    delta->old_hex_index  = (hex_index)-1;
    delta->old_bound_mask = (1 << 4) - 1;

    solver->no_goods_size = 2;
  }
  else {
    delta->old_hex_index  = solver->no_goods[solver->no_goods_size - 2];
    delta->old_bound_mask = solver->no_goods[solver->no_goods_size - 1];
  }

  solver->no_goods_size -= 2;

  uint8_t bound_mask = delta->old_hex_index == hex_id ?
    delta->old_bound_mask : 0;

  if (delta->old_hex_index != hex_id) {
    for (quad_vertex_index i = 0; i < 4; i++) {
      if ((delta->old_bound_mask & (1 << i)) == 0) {
        if (no_goods_reserve(solver, 3) != SUCCESS)
          goto out_of_memory;

        solver->no_goods[solver->no_goods_size++] = i;
        solver->no_goods[solver->no_goods_size++] =
          solver->hexes[8 * delta->old_hex_index + 4 + i];
        solver->no_goods[solver->no_goods_size++] = NO_GOODS_DOWN;
      }
    }

    for (hex_index i = delta->old_hex_index + 1; i <= hex_id; i++) {
      if (no_goods_reserve(solver, 5 + 4 * 3) != SUCCESS)
        goto out_of_memory;

      solver->no_goods[solver->no_goods_size++] =
        NO_GOODS_FACE;
      for (quad_vertex_index j = 0; j < 4; j++)
        solver->no_goods[solver->no_goods_size++] =
          solver->hexes[8 * i + j];

      if (i == hex_id) break;

      for (quad_vertex_index j = 0; j < 4; j++) {
        solver->no_goods[solver->no_goods_size++] = j;
        solver->no_goods[solver->no_goods_size++] =
          solver->hexes[8 * i + 4 + j];
        solver->no_goods[solver->no_goods_size++] = NO_GOODS_DOWN;
      }
    }
  }

  for (quad_vertex_index i = 0; i < 4; i++) {
    if ((bound_mask & (1 << i)) == 0 &&
        (solver->bound_mask & (1 << i)) != 0) {
      if (no_goods_reserve(solver, 3) != SUCCESS)
        goto out_of_memory;

      bound_mask |= (1 << i);

      solver->no_goods[solver->no_goods_size++] = i;
      solver->no_goods[solver->no_goods_size++] =
        solver->hexes[8 * hex_id + 4 + i];
      solver->no_goods[solver->no_goods_size++] = NO_GOODS_DOWN;
    }
  }

  if (no_goods_reserve(solver, 2) != SUCCESS)
    goto out_of_memory;

  solver->no_goods[solver->no_goods_size++] = hex_id;
  solver->no_goods[solver->no_goods_size++] = bound_mask;

  delta->size = solver->no_goods_size - old_size;

  return SUCCESS;

out_of_memory:
  delta->size = 0;
  solver->no_goods_size = old_size;
  solver_undo_no_goods(solver, delta);
  return OUT_OF_MEMORY;
}

error_code solver_is_no_good(solver *solver,
                             quad_vertex_index i, vertex_index v,
                             no_goods_delta *delta) {
  if (solver->num_symmetries == 0)
    return SUCCESS;

  if (solver->no_goods_size == 0)
    return SUCCESS;

  if (no_goods_reserve(solver, 2) != SUCCESS)
    return OUT_OF_MEMORY;

  memmove(solver->no_goods + solver->no_goods_size,
          solver->no_goods + solver->no_goods_size - 2,
          2 * sizeof(*solver->no_goods));

  solver->no_goods[solver->no_goods_size - 2] = i;
  solver->no_goods[solver->no_goods_size - 1] = v;

  solver->no_goods_size += 2;
  delta->size += 2;

  return SUCCESS;
}

error_code solver_hex_is_no_good(solver *solver, no_goods_delta *delta) {
  if (solver->num_symmetries == 0)
    return SUCCESS;

  if (no_goods_reserve(solver, 9) != SUCCESS)
    return OUT_OF_MEMORY;

  solver->no_goods[solver->no_goods_size] = NO_GOODS_HEX;
  memcpy(solver->no_goods + solver->no_goods_size + 1,
         solver->hexes + 8 * solver->hex_id,
         8 * sizeof(*solver->hexes));

  solver->no_goods_size += 9;
  delta->size += 9;

  return SUCCESS;
}

void solver_undo_no_goods(solver *solver, const no_goods_delta *delta) {
  if (solver->num_symmetries != 0) {
    solver->no_goods_size -= delta->size;

    if (solver->no_goods_size != 0 && (solver->flags & SOLVER_SHELLABLE_ONLY) == 0) {
      solver->no_goods[solver->no_goods_size - 2] = delta->old_hex_index;
      solver->no_goods[solver->no_goods_size - 1] = delta->old_bound_mask;
    }
  }
}

/**
 * Symmetry breaking by domination testing. For each symmetry permutation, the
 * algorithm checks if the current state is dominated by one of the no-goods.
 *
 * We use `permutation` to refer to the transformation from the no-good to the
 * current state and `inv_permutation` for the inverse transformation.
 *
 * We use VERTEX_SET_MAX to represent internal vertices whose mapping is not yet
 * known.
 *
 * In order to test for dominance, the algorithm iterates over the hexahedra in
 * the no-good, trying to identify a bijection between the vertices of both
 * solutions.
 *
 * The mapping for the first four vertices in each hex of the no-good is always
 * known, because that face is either on the boundary or has already
 * occurred. This allows the hexahedron containing this face in the other
 * solution to be identified. The mapping between the remaining four vertices is
 * extended if possible (if it isn't, the current state is not dominated by the
 * no-good).
 */

/**
 * Walks through the tree of no-goods to determine if any of them dominate the
 * current node when the given permutation is applied.
 */
static bool is_dominated_by(solver *solver, hex_index hex_id,
                            vertex_index *permutation,
                            vertex_index *inv_permutation);

bool is_dominated(solver *solver) {
  hex_index hex_id = solver->hex_id;
  if (solver->num_symmetries == 0 ||
      solver->no_goods_size == 0) return false;

  vertex_index permutation[VERTEX_SET_MAX];
  vertex_index inv_permutation[VERTEX_SET_MAX];

  for (quad_vertex_index i = 0; i < 4; i++) {
    if ((solver->bound_mask & (1u << i)) == 0)
      solver->hexes[8 * hex_id + 4 + i] = VERTEX_SET_MAX;
  }

  for (size_t i = 0; i < solver->num_symmetries; i++) {
    memcpy(permutation,
           solver->boundary_group + i * solver->num_boundary_vertices,
           solver->num_boundary_vertices * sizeof(*permutation));
    for (vertex_index v_i = solver->num_boundary_vertices;
         v_i < solver->num_vertices; v_i++) {
      permutation[v_i] = VERTEX_SET_MAX;
      inv_permutation[v_i] = VERTEX_SET_MAX;
    }

    invert_permutation(inv_permutation, permutation,
                       solver->num_boundary_vertices);

    if (is_dominated_by(solver, hex_id, permutation, inv_permutation))
      return true;
  }

  return false;
}

void incremental_symmetry_data_trail(const solver *solver, uint8_t *trail) {
  if (solver->num_symmetries == 0) return;

  *(uint64_t*)trail = *(uint64_t*)solver->symmetry_data;

  size_t write_offset = sizeof(uint64_t), read_offset = sizeof(uint64_t);

  size_t entry_size = solver_incremental_symmetry_entry_size(
    solver->num_vertices, solver->num_hexes);

  for (size_t i = 0; i < solver->num_symmetries; i++, read_offset += entry_size) {
    if (*(uint64_t*)trail & (1ul << i)) {
      memcpy(trail + write_offset, solver->symmetry_data + read_offset,
             entry_size);
      write_offset += entry_size;
    }
  }
}

void incremental_symmetry_data_restore(solver *solver, const uint8_t *trail) {
  if (solver->num_symmetries == 0) return;

  *(uint64_t*)solver->symmetry_data = *(uint64_t*)trail;

  size_t write_offset = sizeof(uint64_t), read_offset = sizeof(uint64_t);

  size_t entry_size = solver_incremental_symmetry_entry_size(
    solver->num_vertices, solver->num_hexes);

  for (size_t i = 0; i < solver->num_symmetries; i++, write_offset += entry_size) {
    if (*(uint64_t*)trail & (1ul << i)) {
      memcpy(solver->symmetry_data + write_offset, trail + read_offset,
             entry_size);
        read_offset += entry_size;
    }
  }
}

static bool incremental_is_dominated_by(
  solver *solver, hex_index hex_id, size_t symmetry_id,
  uint64_t *offset, uint64_t *to_check,
  vertex_index *permutation,
  vertex_index *inv_permutation);

bool incremental_domination_test(solver *solver) {
  size_t offset = sizeof(uint64_t);

  for (size_t i = 0; i < solver->num_symmetries; i++) {
    vertex_index *permutation = (vertex_index*)(solver->symmetry_data + offset);
    offset += solver->num_vertices * sizeof(vertex_index);

    vertex_index *inv_permutation = solver->symmetry_data + offset;
    offset += solver->num_vertices * sizeof(vertex_index);

    offset = (offset + alignof(uint64_t) - 1) & ~((size_t)(alignof(uint64_t) - 1));

    uint64_t *index = (uint64_t*)(solver->symmetry_data + offset);
    offset += sizeof(uint64_t);

    uint64_t *to_check = (uint64_t*)(solver->symmetry_data + offset);
    offset += sizeof(uint64_t) * ((solver->num_hexes + 63) / 64);

    if (*(uint64_t*)solver->symmetry_data & (1ul << i)) {
      if (incremental_is_dominated_by(
            solver, solver->hex_id, i,
            index, to_check, permutation, inv_permutation))
        return true;
    }
  }

  return false;
}

static void invert_permutation(vertex_index *inv,
                               const vertex_index *permutation,
                               vertex_index n) {
  for (vertex_index i = 0; i < n; i++)
    inv[permutation[i]] = i;
}

/**
 * Finds the hexahedron which contains a quadrangle, and returns the four
 * vertices opposite to this quadrangle.
 *
 * @param to_check Array containing the indices of the hexahedra to test
 * @param num_check Number of hexahedra to search
 *
 * @param to_check bitset representing the set of already tested hexahedra
 * @param quad The quadrangle to search for
 * @param [out] found_id Index of the hexahedron containing the quadrangle
 * @param [out] opposite The four vertices opposite to the input quadrangle
 *
 * @return True if the quadrangle was found.
 */
static bool find_hex_quad(solver *solver, hex_index hex_id,
                          const uint64_t *to_check,
                          const vertex_index *quad,
                          hex_index *found_id,
                          vertex_index *opposite);

static bool can_reverse_shell_rest(solver *solver, const uint64_t *to_check);

static bool is_dominated_by(solver *solver, hex_index hex_id,
                            vertex_index *permutation,
                            vertex_index *inv_permutation) {
  uint64_t to_check[4] = {0};

  hex_index cur_hex = NO_HEXAHEDRON;
  vertex_index opposite_quad[4];

  size_t offset = 0;
  size_t limit = solver->no_goods_size;
  if ((solver->flags & SOLVER_SHELLABLE_ONLY) == 0) limit -= 2;

  while (offset != limit) {
    uint8_t command = solver->no_goods[offset++];
    if (command < 4) {
      vertex_index in_nogood = solver->no_goods[offset++];
      vertex_index in_state = opposite_quad[command];

      bool followed_by_branch = offset != solver->no_goods_size - 2 &&
        solver->no_goods[offset] >= 4;

      bool success = extend_permutation(permutation, inv_permutation,
                                        in_nogood, in_state);

      /* If we matched a leaf, the current node is dominated. */
      if (success && !followed_by_branch)
        return true;

      /* If we fail to match a branching node, the current node cannot be
       * dominated. */
      if (!success && followed_by_branch)
        return false;
    }
    else if (command == NO_GOODS_FACE) {
      /* Find the hexahedra that contains (the permuted form of) the first
       * face of the next hexahedron. */
      vertex_index in_quad[4];
      for (quad_vertex_index i = 0; i < 4; i++)
        in_quad[i] = permutation[solver->no_goods[offset++]];

      if (!find_hex_quad(solver, hex_id, to_check,
                         in_quad, &cur_hex,
                         opposite_quad))
        return false;

      to_check[cur_hex / 64] |= (1ul << (cur_hex % 64));
    }
    else if (command == NO_GOODS_HEX) {
      bool followed_by_branch = offset + 8 != solver->no_goods_size &&
        solver->no_goods[offset + 8] != NO_GOODS_HEX;

      bool can_reverse_shell = true;

      vertex_index in_quad[4];
      for (quad_vertex_index i = 0; i < 4; i++)
        in_quad[i] = permutation[solver->no_goods[offset++]];

      bool success = true;

      if (!find_hex_quad(solver, hex_id, to_check,
                         in_quad, &cur_hex,
                         opposite_quad)) {
        success = false;
      }
      else {
        to_check[cur_hex / 64] |= (1ul << (cur_hex % 64));

        uint8_t new_mask = 0;
        quad_vertex_index i;
        for (i = 0; i < 4; i++) {
          vertex_index in_state  = opposite_quad[i];
          vertex_index in_nogood = solver->no_goods[offset + i];

          bool is_new = inv_permutation[in_state] == VERTEX_SET_MAX;

          if (!extend_permutation(permutation, inv_permutation,
                                  in_nogood, in_state))
            break;

          if (is_new) new_mask |= (1 << i);
        }

        if (i == 4 && !followed_by_branch)
          can_reverse_shell = can_reverse_shell_rest(solver, to_check);

        if (i != 4 || !can_reverse_shell) {
          success = false;

          for (quad_vertex_index j = 0; j < i; j++) {
            if (new_mask & (1 << j)) {
              permutation[solver->no_goods[offset + j]] = VERTEX_SET_MAX;
              inv_permutation[opposite_quad[j]] = VERTEX_SET_MAX;
            }
          }

          to_check[cur_hex / 64] &= ~(1ul << (cur_hex % 64));
        }
      }

      offset += 4;

      if (success && !followed_by_branch)
        return true;

      if (!success && followed_by_branch)
        return false;
    }
  }

  return false;
}

static bool find_hex_quad(solver *solver, hex_index hex_id,
                          const uint64_t *to_check,
                          const vertex_index *quad,
                          hex_index *found_id,
                          vertex_index *opposite) {
  quad_descriptor desc = solver->quad_descriptors[
    diagonal_index(solver->num_vertices, quad[0], quad[2])];

  quad_index quad_id = quad_descriptor_quad_index(desc);

  if (solver->num_occurrences[quad_id] != 0)
    return false;

  for (hex_index i = 0; i <= hex_id; i++) {
    if (to_check[i / 64] & (1ul << (i % 64)))
      continue;

    vertex_index *hex = solver->hexes + 8 * i;
    hex_vertex_index loc[4];

    quad_vertex_index j;
    for (j = 0; j < 4; j++) {
      for (loc[j] = 0; loc[j] < 8 && hex[loc[j]] != quad[j]; loc[j]++);
      if (loc[j] == 8) break;
    }

    if (j != 4) continue;

    *found_id = i;

    for (j = 0; j < 4; j++) {
      hex_vertex_index n_a = HEX_NEIGHBORS[loc[j]][0];
      hex_vertex_index n_b = HEX_NEIGHBORS[loc[j]][1];
      hex_vertex_index n_c = HEX_NEIGHBORS[loc[j]][2];

      if (n_a != loc[0] && n_a != loc[1] && n_a != loc[2] && n_a != loc[3])
        opposite[j] = hex[n_a];
      else if (n_b != loc[0] && n_b != loc[1] && n_b != loc[2] && n_b != loc[3])
        opposite[j] = hex[n_b];
      else
        opposite[j] = hex[n_c];
    }

    return true;
  }

  return false;
}

static bool can_create_quad(solver *solver, const vertex_index *quad);

static bool incremental_is_dominated_by(
  solver *solver, hex_index hex_id, size_t symmetry_id,
  uint64_t *offset, uint64_t *to_check,
  vertex_index *permutation,
  vertex_index *inv_permutation) {
  size_t limit = solver->no_goods_size;
  if ((solver->flags & SOLVER_SHELLABLE_ONLY) == 0) limit -= 2;

  vertex_index temp_permutation[solver->num_vertices];
  vertex_index temp_inv_permutation[solver->num_vertices];
  uint64_t temp_offset;
  uint64_t temp_to_check[4];

  bool on_temp_data = false;

  while (*offset != limit) {
    uint8_t command = solver->no_goods[(*offset)++];
    if (command == NO_GOODS_HEX) {
      bool followed_by_branch = (*offset) + 8 != solver->no_goods_size &&
        solver->no_goods[(*offset) + 8] != NO_GOODS_HEX;

      bool can_reverse_shell = true;

      vertex_index in_quad[4];
      for (quad_vertex_index i = 0; i < 4; i++)
        in_quad[i] = permutation[solver->no_goods[(*offset)++]];

      bool success = true;

      hex_index cur_hex = NO_HEXAHEDRON;
      vertex_index opposite_quad[4];
      if (!find_hex_quad(solver, hex_id, to_check,
                         in_quad, &cur_hex,
                         opposite_quad)) {
        if (followed_by_branch) {
          *offset -= 5;
          return false;
        }
        else {
          if (!on_temp_data && !can_create_quad(solver, in_quad)) {
            on_temp_data = true;

            memcpy(temp_permutation, permutation, sizeof(temp_permutation));
            memcpy(temp_inv_permutation, inv_permutation, sizeof(temp_inv_permutation));
            memcpy(temp_to_check, to_check, sizeof(uint64_t) * ((solver->num_hexes + 63) / 64));
            temp_offset = *offset;

            *offset -= 5;

            permutation = temp_permutation;
            inv_permutation = temp_inv_permutation;
            to_check = temp_to_check;
            offset = &temp_offset;
          }

          success = false;
        }
      }
      else {
        to_check[cur_hex / 64] |= (1ul << (cur_hex % 64));

        uint8_t new_mask = 0;
        quad_vertex_index i;
        for (i = 0; i < 4; i++) {
          vertex_index in_state  = opposite_quad[i];
          vertex_index in_nogood = solver->no_goods[(*offset) + i];

          bool is_new = inv_permutation[in_state] == VERTEX_SET_MAX;

          if (!extend_permutation(permutation, inv_permutation,
                                  in_nogood, in_state))
            break;

          if (is_new) new_mask |= (1 << i);
        }

        if (i == 4 && !followed_by_branch)
          can_reverse_shell = can_reverse_shell_rest(solver, to_check);

        if (i != 4 || !can_reverse_shell) {
          success = false;

          for (quad_vertex_index j = 0; j < i; j++) {
            if (new_mask & (1 << j)) {
              permutation[solver->no_goods[(*offset) + j]] = VERTEX_SET_MAX;
              inv_permutation[opposite_quad[j]] = VERTEX_SET_MAX;
            }
          }

          to_check[cur_hex / 64] &= ~(1ul << (cur_hex % 64));
        }
      }

      *offset += 4;

      if (success && !followed_by_branch) {
        *offset -= 9;
        return true;
      }

      if (!success && followed_by_branch) {
        if (!on_temp_data) {
          *(uint64_t*)solver->symmetry_data &= ~(1ul << symmetry_id);
          *offset -= 9;
        }
        return false;
      }
    }
  }

  return false;
}

static bool can_create_quad(solver *solver, const vertex_index *quad) {
  if (solver->quad_descriptors[diagonal_index(
          solver->num_vertices, quad[0], quad[2])] != EMPTY_QUAD_DESCRIPTOR)
    return false;

  if (solver->quad_descriptors[diagonal_index(
          solver->num_vertices, quad[1], quad[3])] != EMPTY_QUAD_DESCRIPTOR)
    return false;

  for (quad_vertex_index i = 0; i < 4; i++) {
    if (vertex_set_contains(&solver->finished_vertices, quad[i]))
      return false;
  }

  return true;
}

static
hex_index initial_mesh_data(solver *solver, const uint64_t *to_check,
                            vertex_set *edges,
                            vertex_set *quad_diagonals) {
  for (vertex_index i = 0; i < solver->max_vertex_id; i++) {
    vertex_set_init_empty(edges + i);
    vertex_set_init_empty(quad_diagonals + i);
  }

  uint64_t on_boundary[4] = {0};

  for (uint32_t i = 0; i < solver->num_boundary_vertices - 2; i++)
    on_boundary[i / 64] |= (1ull << (i % 64));

  hex_index num_seen = 0;
  for (hex_index i = 0; i <= solver->hex_id; i++) {
    if (to_check[i / 64] & (1ul << (i % 64))) {
      vertex_index *hex = solver->hexes + 8 * i;

      for (hex_facet_index j = 0; j < 6; j++) {
        vertex_index a = hex[HEX_FACET_CONTENTS[j][0]];
        vertex_index b = hex[HEX_FACET_CONTENTS[j][2]];

        quad_descriptor desc = solver->quad_descriptors[
          diagonal_index(solver->num_vertices, a, b)];
        quad_index q = quad_descriptor_quad_index(desc);

        on_boundary[q / 64] ^= (1ull << (q % 64));
      }
      num_seen++;
    }
  }

  for (quad_index i = 0; i < solver->max_quad_id; i++) {
    if ((on_boundary[i / 64] & (1ull << (i % 64))) == 0) continue;
    const vertex_index *quad = solver->quads + 4 * i;
    for (quad_vertex_index j = 0; j < 4; j++) {
      vertex_set_insert(&edges[quad[j]], quad[(j + 1) % 4]);
      vertex_set_insert(&edges[quad[(j + 1) % 4]], quad[j]);

      vertex_set_insert(&quad_diagonals[quad[j]], quad[(j + 2) % 4]);
    }
  }

  return num_seen;
}

typedef enum intersection_mask_bits {
  INTERSECTION_MASK_BIT_VERTEX0 = 0,
  INTERSECTION_MASK_BIT_EDGE0   = 8,
  INTERSECTION_MASK_BIT_FACET0  = 20,
} intersection_mask_bits;

static const uint32_t FACET_INTERSECTION_MASKS[] = {
  0x0102b0f, /* 0b000001_000000101011_00001111 */
  0x0221533, /* 0b000010_001000010101_00110011 */
  0x0445866, /* 0b000100_010001011000_01100110 */
  0x088e0cc, /* 0b001000_100011100000_11001100 */
  0x1018699, /* 0b010000_000110000110_10011001 */
  0x20f00f0, /* 0b100000_111100000000_11110000 */
};

static const uint32_t FACE_SELECTOR = ((1u << 6)-1) << INTERSECTION_MASK_BIT_FACET0;

static bool can_peel_hex(
  const vertex_set *edges, const vertex_set *quad_diagonals,
  const vertex_index *hex, uint32_t *intersection_mask) {
  uint32_t mask = 0, expected_mask = 0;
  for (hex_vertex_index i = 0; i < 8; i++) {
    if (!vertex_set_empty(&edges[hex[i]]))
      mask |= (1ul << (INTERSECTION_MASK_BIT_VERTEX0 + i));
  }

  for (hex_edge_index i = 0; i < 12; i++) {
    vertex_index a = hex[HEX_EDGES[i][0]];
    vertex_index b = hex[HEX_EDGES[i][1]];
    if (vertex_set_contains(&edges[a], b))
      mask |= (1ul << (INTERSECTION_MASK_BIT_EDGE0 + i));
  }

  uint32_t num_faces = 0;

  for (hex_facet_index i = 0; i < 6; i++) {
    vertex_index a = hex[HEX_FACET_CONTENTS[i][0]];
    vertex_index b = hex[HEX_FACET_CONTENTS[i][2]];

    if (vertex_set_contains(&quad_diagonals[a], b)) {
      mask |= (1ul << (INTERSECTION_MASK_BIT_FACET0 + i));
      expected_mask |= FACET_INTERSECTION_MASKS[i];
      num_faces++;
    }
  }

  *intersection_mask = mask;

  if (mask == 0 || mask != expected_mask)
    return false;

  if (num_faces == 2 || num_faces == 4) {
    uint32_t face_mask = (mask & FACE_SELECTOR) >> INTERSECTION_MASK_BIT_FACET0;
    if (num_faces == 4) face_mask = ~face_mask;

    hex_facet_index f_a = find_lsb32(face_mask);
    face_mask &= ~(1u << f_a);

    hex_facet_index f_b = find_lsb32(face_mask);

    return (FACET_INTERSECTION_MASKS[f_a] & FACET_INTERSECTION_MASKS[f_b]) != 0;
  }
  else
    return true;
}

static void peel_hex(vertex_set *edges, vertex_set *quad_diagonals,
                     const vertex_index *hex, uint32_t intersection_mask) {
  for (hex_edge_index i = 0; i < 12; i++) {
    if ((intersection_mask & (1 << (INTERSECTION_MASK_BIT_EDGE0 + i))) == 0) {
      vertex_index a = hex[HEX_EDGES[i][0]];
      vertex_index b = hex[HEX_EDGES[i][1]];

      vertex_set_insert(&edges[a], b);
      vertex_set_insert(&edges[b], a);
    }
  }

  for (hex_facet_index i = 0; i < 6; i++) {
    if ((intersection_mask & (1 << (INTERSECTION_MASK_BIT_FACET0 + i))) == 0) {
      for (quad_vertex_index j = 0; j < 4; j++) {
        vertex_index a = hex[HEX_FACET_CONTENTS[i][j]];
        vertex_index b = hex[HEX_FACET_CONTENTS[i][(j + 2) % 4]];

        vertex_set_insert(&quad_diagonals[a], b);
      }
    }
  }
}

static void reinsert_hex(vertex_set *edges, vertex_set *quad_diagonals,
                         const vertex_index *hex, uint32_t intersection_mask) {
  for (hex_edge_index i = 0; i < 12; i++) {
    if ((intersection_mask & (1 << (INTERSECTION_MASK_BIT_EDGE0 + i))) == 0) {
      vertex_index a = hex[HEX_EDGES[i][0]];
      vertex_index b = hex[HEX_EDGES[i][1]];

      vertex_set_remove(&edges[a], b);
      vertex_set_remove(&edges[b], a);
    }
  }

  for (hex_facet_index i = 0; i < 6; i++) {
    if ((intersection_mask & (1 << (INTERSECTION_MASK_BIT_FACET0 + i))) == 0) {
      for (quad_vertex_index j = 0; j < 4; j++) {
        vertex_index a = hex[HEX_FACET_CONTENTS[i][j]];
        vertex_index b = hex[HEX_FACET_CONTENTS[i][(j + 2) % 4]];

        vertex_set_remove(&quad_diagonals[a], b);
      }
    }
  }
}

#define MAX_REVERSE_SHELL_DISCREPENCY 10

static bool try_reverse_shell(
  solver *solver, uint64_t *to_check,
  vertex_set *edges, vertex_set *quad_diagonals,
  hex_index num_seen, size_t *discrepency) {
  if (num_seen == solver->hex_id + 1)
    return true;

  for (hex_index i = 0; i <= solver->hex_id && (*discrepency) <= MAX_REVERSE_SHELL_DISCREPENCY; i++) {
    if (to_check[i / 64] & (1ul << (i % 64)))
      continue;

    uint32_t intersection_mask;
    if (can_peel_hex(edges, quad_diagonals,
                     solver->hexes + 8 * i, &intersection_mask)) {
      peel_hex(edges, quad_diagonals, solver->hexes + 8 * i, intersection_mask);
      to_check[i / 64] |= (1ul << (i % 64));

      if (try_reverse_shell(solver, to_check, edges, quad_diagonals, num_seen + 1,
                            discrepency))
        return true;
      else
        (*discrepency)++;

      to_check[i / 64] &= ~(1ul << (i % 64));
      reinsert_hex(edges, quad_diagonals, solver->hexes + 8 * i, intersection_mask);
    }
  }

  return false;
}

static bool can_reverse_shell_rest(solver *solver, const uint64_t *to_check) {
  vertex_set edges[VERTEX_SET_MAX];
  vertex_set quad_diagonals[VERTEX_SET_MAX];
  uint64_t bitset[4];
  memcpy(bitset, to_check, sizeof(*bitset) * (solver->num_hexes + 63) / 64);

  hex_index num_seen = initial_mesh_data(solver, to_check,
                                         edges, quad_diagonals);

  size_t discrepency = 0;

  return try_reverse_shell(
    solver, bitset,
    edges, quad_diagonals, num_seen, &discrepency);
}
